Quantum confined electronic states in atomically well-defined graphene nanostructures 
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Despite the enormous interest in the properties of graphene and the potential of graphene nanos- 
tructures in electronic applications, the study of quantum confined states in atomically well-defined 
graphene nanostructures remains an experimental challenge. Here, we study graphene quantum 
dots (GQDs) with well-defined edges in the zigzag direction, grown by chemical vapor deposition 
(CVD) on an iridium(lll) substrate, by low-temperature scanning tunneling microscopy (STM) 
and spectroscopy (STS). We measure the atomic structure and local density of states (LDOS) of 
individual GQDs as a function of their size and shape in the range from a couple of nanometers 
up to ca. 20 nm. The results can be quantitatively modeled by a relativist ic wave equation and 
atomistic tight-binding calculations. The observed states are analogous to the solutions of the text 
book "particle-in-a-box" problem applied to relativistic massless fermions. 



PACS numbers: 73.21. La 73.22.Pr 73.63.Kv 81.05.ue 

Graphene, a monolayer of carbon atoms that is a 2-D 
metal where the charge carriers behave as massless rel- 
ativistic electrons has attracted enormous scientific and 
technological interest. [1-4]. Despite the potential of 
graphene nanostructures in electronic applications [5- 
12], the study of quantum confined electronic states in 
atomically well-defined graphene nanostructures remains 
an experimental challenge. Basic questions, such as the 
relation between the atomic configuration of graphene 
nanostructures and the spatial distribution and energy 
of their electronic states have not been experimentally 
addressed. 

In previous experiments, macroscopic graphene sheets 
have been studied by scanning tunneling microscopy 
(STM) and spectroscopy (STS), focusing on the elec- 
tronic structure and scattering processes in epitaxial 
graphene [13, 14] and the density of states and charge 
puddles in graphene sheets deposited on an insulator 
[15-18]. It is, however, also possible to grow much 
smaller graphene nanostructures (graphene quantum 
dots, GQDs) by chemical vapor deposition (CVD), and 
characterize them with scanning probe methods [10, 19- 
21]- 

In this report, we present low-temperature STM and 
STS experiments on GQDs with well-defined atomic 
structures grown by CVD on an Ir(lll) substrate. We 
can readily access individual GQDs and measure their 
atomic structure with STM. Measurement of the local 
density of states (LDOS, proportional to the dl/dV^ sig- 
nal) allows us to probe the spatial structure and energy 
of the quantum confined energy levels for GQDs with 
variable size and shape. The measured LDOS maps can 
be reproduced by tight-binding (TB) calculations, where 
we use the exact atomic structure of the GQDs as deter- 



mined by STM, and by the Klein-Gordon (KG) equation, 
which is a continuum model describing particles with lin- 
ear dispersion. 

The Ir(lll) surface was cleaned by sputtering at 1100 
K and annealing at 1500 K. After the sample had cooled 
below 570 K, ethylene was deposited (3 x 10 -6 mbar for 
10 s) on the surface. The GQD size could then be con- 
trolled by the growth temperature [21]: larger (smaller) 
GQDs were grown by heating the sample to 1470 K (1170 
K) for 10 s. After the CVD growth of the GQDs, the 
sample was inserted into a low-temperature STM (T = 
4.8 K, Omicron LT-STM) housed within the same ultra- 
high vacuum system (base pressure < 10 -10 mbar). We 
used cut Ptlr tips and the bias voltage (Vb) was defined 
as sample voltage with respect to the tip. The dl/dV^ 
signal was recorded with a lock-in amplifier by applying 
a small sinusoidal variation to the bias voltage (typically 
30 mV rms at 660 Hz). This gives an energy resolution 
of ca. 75 meV in our experiments [22]. To ensure that 
the modulation would not couple to the feedback loop, 
a 300 Hz low-pass filter was used for the feedback in- 
put. The experimental dl/dVh images are averages of 
the trace and retrace scans (Fig. 1) and of the trace and 
retrace scans of two consecutive images (up and down, 
Fig. 2) to increase the signal to noise ratio. 

Fig. 1(a) shows a large-scale overview scan of a typical 
sample. We find interconnected graphene patches (indi- 
cated by G) as well as small isolated GQDs (red dotted 
circles). The CVD growth yields a relatively broad dis- 
tribution of different GQD sizes ranging from a couple 
of nanometers up to ca. 20 nm, most of them with a 
roughly hexagonal shape. All the GQDs have edges in 
the zigzag direction (corresponding to the close-packed 
atomic rows of the underlying Ir(lll) surface) with a 
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FIG. 1. (Color online) STM imaging and spectroscopy on 
GQDs grown by CVD on Ir(lll). (a) Large scale STM im- 
age of graphene islands (G) on an iridium(lll) substrate (ac- 
quired at / = 40 pA and Vb = 1.0 V). Small graphene QDs 
have been indicated by red circles, (b) STM topography of a 
small GQD (-0.05 V / 100 pA) with an overlaid atomic model 
which has perfect hexagonal symmetry with 7 benzene ring 
long edges, (c) dln//dlnVb spectra measured on the points 
indicated in (b), the green lines indicate the bias voltages 
corresponding to the LDOS maps shown in panel (d). (d) 
Measured LDOS maps (gray line denotes the edges of the 
GQD) at bias voltages corresponding to the two resonances 
in the spectrum shown in panel (c). (e) The corresponding 
LDOS maps calculated for a particle in a box at the indicated 
energies and the underlying eigenstates. 

very small roughness. We see kinks of one or two atomic 
rows at the GQD edges [23]. Closer inspection of small 
GQDs at a bias voltage close to zero bias shows that 
the edges are bright both in the actual STM topogra- 
phy as well as in the simultaneously recorded dl/dV^ 
images. These edge states are expected for zigzag edges 
in graphene [2, 24-26]. More information can be found 
in the Supplementary Information [23] . 

We now focus on the delocalized, quantum confined 
states inside the GQDs. We can map the atomic struc- 
ture of the GQD by STM as shown for a small GQD 
with perfect hexagonal symmetry with 7 benzene ring 
long edges in Fig. 1(b). The LDOS can be accessed 
through dlnl /dlnVh measurements as shown in Fig. 
1(c); we clearly observe an increased and spatially de- 
pendent LDOS on the GQDs. There is a pronounced 
maximum of the LDOS measured in the center of the 
GQD (blue line in Fig. 1(c)) at a bias of -0.6 V. Mov- 
ing away from the center of the GQD, the intensity of 
this peak is reduced, and another resonance emerges at 



a bias of -0.9 V (red line). We can map the spatial 
shape of the orbitals responsible for these resonances by 
measuring the dl/dVh signal during STM imaging under 
constant-current feedback at biases corresponding to the 
resonances [Fig. 1(d)]. These states have the familiar 
appearance of the lowest energy levels of the text book 
particle-in-a-box problem and can be characterized us- 
ing symmetry labels borrowed from atomic physics. The 
lowest state has IS symmetry (no nodal planes) and the 
first excited state is composed of two IP type orbitals 
(1P X and lP y ) which are degenerate in this case of a per- 
fect hexagonal GQD. STM probes the sum of the squared 
wavefunctions ip 2 P + ip 2 P leading to a doughnut shaped 
dl/dVh signal as we observe in the experiment. Compar- 
ison of these states with TB calculations can be found in 
the Supplementary Information [23]. 

We note here that at positive bias, electronic reso- 
nances with clear peaks in the dl/dV^ spectrum cannot 
be observed [23]. Based on DFT calculations on Ir(lll), 
there is a dense set of energy bands above the Fermi 
energy at the K point of Brillouin zone. It is likely that 
interaction with these states masks the intrinsic graphene 
states at positive bias [27]. 

These experiments can be reproduced by both TB cal- 
culations and by a continuum model for particles with 
linear dispersion confined to a GQD [23]. Here we use 
the KG equation [28, 29] 

- v 2 n 2 v 2 ^ = Efyi (i) 

where is the Fermi velocity (10 6 m/s in isolated 
graphene) and the boundary condition is given by ipi = 
at the edges of the GQD. A more accurate boundary con- 
dition would be needed to take into account the sublattice 
pseudospin and the interaction with the Ir substrate. It 
is clear that the KG equation cannot be used to model 
the edge states (in contrast to the Dirac equation and TB 
calculations [2, 24, 25]). However, as shown below, the 
LDOS plots from the KG equation are remarkably simi- 
lar to the TB calculations and the experimental results, 
although the number of states in a given energy interval 
is too small. We use the experimentally determined ge- 
ometries of the GQD in our calculations [23]. The lowest 
energy solutions of Eq. (1) are plotted in Fig. 1(e) as the 
squared wavefunctions corresponding to the experimen- 
tally measured LDOS ex ^2 SE i/j 2 , where SE is the energy 
resolution of the experiment [30] . 

We have measured the LDOS at different bias voltages 
on a larger GQD shown in Fig. 2(a). The periodic vari- 
ation with a period of 2.5 nm seen on the topographic 
STM images is a moire pattern resulting from the lattice 
mismatch between graphene and Ir [14, 21]. The STM 
contrast results mostly from a small (ca. 30 pm) geo- 
metric modulation of the graphene structure [14]. Our 
calculations neglecting this moire-induced potential mod- 
ulation yield quantitative agreement with the experiment 
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FIG. 2. (Color online) Detailed comparison between STM and 
STS experiments and computational results on a large GQD. 
(a) Atomically resolved STM image of the GQD (7 = 3 nA, Vb 
= 1 mV). (b) d//dVb maps recorded under constant-current 
STM feedback at the bias voltages indicated in the figure (I 
= 1 nA). (c,d) Corresponding LDOS plots at the indicated 
energies calculated using a TB model (c) and the KG equation 
(d) as described in the text, (e) Correspondence between the 
experimental and the calculated energies based on TB (red 
squares) and the KG equation (blue circles) calculated with 
v F 6.2 x 10 5 m/s. 



and the expected potential modulation due to the moire 
pattern is small compared to the confinement energy in 
our GQDs. It has been reported that the size and shape 
of the GQDs is influenced by the moire pattern and the 
edges prefer to run along the fee and hep regions of the 
moire [21, 31]. We also observe GQDs that are smaller 
than the moire period (6x6 and 7x7). For larger GQDs, 
the kinks on the edges are spaced by roughly one moire 
period. 

The asymmetry of the GQD breaks the degeneracies 
(e.g. 1P X and lP y states) of the purely hexagonal GQD. 
This can be seen in the measured LDOS maps shown in 
Fig. 2(b) (The Ir substrate has been removed in the im- 
ages using the simultaneously acquired STM topography 
image as a mask, images with the background can be 
found in the Supplementary Information [23]): after the 
IS state (bias -0.25 V), we observe increased intensity 
at the top and bottom end of the GQD consistent with 
the lP y envelope wavefunction along the long GQD axis 
(at -0.30 V). At more negative bias, the 1P X state also 
contributes and the long GQD edges are brighter (-0.35 
V). Subsequently, the next eigenstate becomes relevant, 
which is seen as an increased intensity in the middle of 



the QD (bias -0.4 V). 

In order to compare experiment and theory in detail, 
we have generated a series of theoretical LDOS maps, 
which are calculated as a weighted and broadened sum 
of squares of TB molecular orbitals (MOs) or KG eigen- 
states close to a given energy [see Figures 2(c,d)] [23]. 
This broadening is justified due to the intrinsic resolution 
of the measurement (75meV) and the life-time broaden- 
ing of the states. In the case of the calculations based on 
the KG equation, the eigenfunctions are given by the so- 
lution of Eq. (1) using the overall shape of the GQD. In 
the TB calculations (we use third-nearest-neighbor TB) 
[2, 24, 32], they correspond to the calculated MOs for 
the GQD with an exact atomic structure as obtained 
from experiment [Fig. 2(a)] [23]. It can be seen that the 
eigenstates of the KG equation (overall geometry) match 
with clusters of TB MOs (exact atomic lattice). Fur- 
thermore, there is a remarkable agreement in how both 
calculated LDOS maps evolve with energy and how the 
experimental conductance maps evolve with the bias. 

Based on a comparison between the experimental and 
computed LDOS maps, we have identified energy / bias 
voltage pairs that give the same spatial features in the 
LDOS with an associated error estimate indicated by er- 
ror bars in Fig. 2(e) [23]. It is clear that with the Fermi- 
velocity f p as the only adjustable parameter (in the case 
of TB calculations, is directly related to the value 
of the hopping integrals), both calculations agree strik- 
ingly well with the experiments. This is also evident from 
Fig. 2(e), where we show the correspondence between the 
experimental bias voltages and the theoretical energies. 
This gives the Fermi velocity = (6.2 ±0. 1) x 10 5 m/s as 
the best-fit to both the KG equation and the TB calcula- 
tions. The two theories yield slightly different values for 
the doping of the GQD, i.e. the intercept of the ?/-axis, 
due to the differences in the theoretical approaches. 

Do we see the peculiar nature of the charge carriers in 
graphene in these LDOS maps? In fact, the Schrodinger 
equation predicts wavefunctions with an identical spatial 
shape as the KG equation since both are second order dif- 
ferential equations; the corresponding eigenenergies are 
related as E$ = E^ G /2mVp. This also explains the dif- 
ferent dispersion relations for free electrons, which are 
either parabolic (Schrodinger) or linear (Klein-Gordon). 
Moreover, the energy of the lowest (and the other) quan- 
tum confined state scales as 1/A 1 / 2 (A is the area of the 
QGD) in the case of the relativistic massless particles, 
instead of 1/A for the particles obeying the Schrodinger 
equation. We demonstrate in Fig. 3 that the charge car- 
riers in our GQDs fulfil the conditions of E oc 1/A 1 / 2 
and have a linear dispersion. Fig. 3(a) shows the bias 
voltage corresponding to the lowest quantum confined 
energy level (determined by the peak position in dl/dVh 
vs. Vb spectra acquired at the center of the GQD) on 
many different GQDs [topographies shown in Fig. 3(b)] 
as a function of the experimentally determined area. The 
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FIG. 3. (Color online) Electronic structure of GQDs as a 
function of their size, (a) STM sample bias corresponding 
to the S state as a function of the area A of the GQD. The 
solid and dashed lines are fits to 1/A 1 ^ 2 and 1/A scaling, 
respectively, (b) Composition of the STM topographies of 
the GQDs used in panel (a) (with different scan parameters), 
(c) A plot of the bias voltages from the STM experiments (x- 
axis) and the energies that give comparable LDOS calculated 
from the Klein-Gordon equation using a single value for vf = 
6.2 x 10 5 m/s (y-axis). 

solid line showing the expected 1/A 1 / 2 scaling fits the 
data clearly better than the 1/A (dashed line) behavior. 

In Fig. 3(c), we present the correspondence between 
experimental bias voltages (x-axis) and the theoretical 
energies calculated with the KG equation (y-axis) for 
many states on several GQDs. The one-to-one correspon- 
dence confirms that the experimental data is consistent 
with the linear dispersion of the Klein-Gordon equation. 
The corresponding Fermi velocity vp = (6.2 ± 0.3) x 10 5 
m/s is slightly smaller than the previous results on 
macroscopic graphene samples on Ir(lll) obtained by 
ARPES (6.5 x 10 5 to 9.2 x 10 5 m/s) [33-35]. Possible 
reasons for this discrepancy are that our STM measure- 
ments probe the average Fermi velocity around the Dirac 
cone (in contrast to ARPES) and our experiments are 
carried out on GQDs instead of bulk graphene. Remark- 
ably, vp remains constant down to the smallest structures 
that we have measured. The intercept with the y-axis in 
Fig. 3(c) and the extrapolation to infinite GQD area in 
Fig. 3(a) indicate that GQDs on Ir(lll) are n-doped by 
- 0.1 eV. 

In summary, we have presented low-temperature STM 
and STS experiments aimed at understanding the quan- 
tum confined energy levels and their spatially resolved 
wavefunctions in atomically well-defined graphene quan- 
tum dots. The measured resonances and corresponding 
LDOS maps correspond to a number of molecular orbitals 
close in energy, calculated by TB for the exact atomic 
geometry. The energy position and LDOS structure of 
these clustered states can also be calculated from the rel- 



ativistic wave-equation for massless particles. Our results 
provide experimental verification of the physics relevant 
for graphene-based opto-electronics where wavefunction 
engineering via well-defined nanostructuring is likely to 
be a central issue. In addition, our experiments indi- 
cate that the intrinsic electronic states of graphene can 
be studied on weakly interacting metal substrates (e.g. 
Ir(lll)). These systems can act as future test beds for 
studying the effects of chemical modifications or doping 
of graphene. 
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Note added -During the review of this letter, we be- 
came aware of related experiments presented in Refs. 
[36, 37]. 
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ORIENTATION AND EDGE STRUCTURE OF GRAPHENE QUANTUM DOTS 




FIG. 1. STM images and their Fourier transforms showing GQD edge orientation and the alignment 
of the Ir and graphene lattices, (a) Atomically resolved STM image of the GQD presented in Fig. 
2 of the main text (5 mV / InA). The contrast is enhanced by an adsorbed molecule on the tip, 
which also inverts the atomic contrast, (b) Square root of the modulus of the FFT of the panel 
(a), (c) Atomic resolution STM image (50 mV / 200 pA) of a larger graphene flake edge where the 
kinks of one atomic row are clearly visible. The inset shows the FFT spectrum of the image. 

The lattices of all the observed graphene islands were roughly aligned with the underlying 
Ir lattice, which is the predominant growth phase of graphene on Ir [1]. This can be clearly 
seen in atomically resolved STM images and their Fourier transforms [Supplementary Figs. 
1(a) and 1(b)]. The Fourier transforms show intensity maxima in a hexagonal pattern 
around the maxima produced by the graphene lattice. These spots arise from the lattice 
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mismatch between graphene and iridium and are analogous to the diffraction pattern of 
graphene on Ir observed in low-energy electron diffraction (LEED) [2]. In our case, these 
hexagons are aligned with the maxima from the graphene lattice indicating that the two 
lattices are aligned. The mismatch in lattice constants is also observed in the real space 
STM images as a moire pattern on graphene islands with a period of around 2.5 nm [2, 3]. 

The orientation of the GQD edges was determined from atomic resolution STM images 
and FFT spectra. All edges were terminated in the zigzag direction. This could be verified 
from FFT spectra of the GQDs, where the first order maxima were always oriented perpen- 
dicular to the graphene flake edges [Supplementary Fig. 1(c)]. The edges always had 120 
degree corners; no 60 degree corners were observed. Thus all GQDs had an overall hexagonal 
shape and no triangular QDs were observed. 

EDGE STATES IN SMALL GRAPHENE QDS 




FIG. 2. STM measurements on edge states in small GQDs. (a) STM topography of an isolated 
GQD flake (/ = 100 pA / V\> = 0.05 V). (b) Simultaneously acquired dl/dVh map under constant- 
current STM feedback. 

Closer inspection of a small GQD [Supplementary Figs. 2(a) and (b)] at a bias voltage 
close to zero bias shows that the edges are bright both in the actual STM topography image 
[Supplementary Fig. 2(a)] and in the simultaneously recorded dl/dV^ image [Supplementary 
Fig. 2(b)]. The enhanced conductance is due to edge states that are expected for zigzag or 
reconstructed zigzag edges [4, 5]. We observe a higher LDOS at the corner and kink sites 
compared to the other edge sites, while atomistic tight-binding modeling predicts highest 
intensity at the middle of the edges with vanishing intensity in the corners of the GQD 
[6]. This disagreement is likely to be related to the remaining coupling of the GQDs with 
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the Ir(lll) substrate. Angle-resolved photoemission (ARPES) experiments have shown that 
close to the Fermi-level, there is an onset of a Ir(lll) surface state that interacts with the 
graphene layer [7, 8]. At negative bias voltages, but before the onset of the delocalized states 
in the interior of the GQD, the experimental results on the edge-state LDOS agree with the 
predictions from tight-binding calculations. Bias-dependent imaging indicates that the edge 
states have a very flat dispersion, i.e. their spatial (exponential) decay away from the edges 
of the GQD is roughly energy independent with a decay constant of ca. 0.5 nm. Generally, 
we can observe the intrinsic properties of GQDs with little interference from the substrate 
in the bias region below -0.1 V. 



COMPARISON OF THE EXPERIMENTS WITH TIGHT-BINDING CALCULA- 
TIONS FOR A SMALL GQD 
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FIG. 3. (a) LDOS calculated by TB compared with experimental results reproduced from Fig. lc 
of the main text. The TB energy levels are indicated by lines at the top of the graph (some levels 
are degenerate). The eigenenergies from the Klein-Gordon equation are indicated by magenta 
lines on the bottom of the graph, (b) Atomic model of the GQD indicating the positions where 
the LDOS is calculated in panel (a), (c) The calculated LDOS map based on TB at the energy 
corresponding to the IS state, (d) calculated LDOS map based on TB at the energy corresponding 
to the IP state. 
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DETERMINATION OF THE GQD SHAPE AND ATOMIC STRUCTURE 





FIG. 4. Examples of the atomic models created for two different GQDs. (a) Atomic model created 
directly from the STM image for a GQD with 6 carbon rings per edge. The grey hexagon around 
the atomic model (middle column) represents the area used for the FEM calculation of the Klein- 
Gordon equation, (b) Example of a larger GQD where the atomic model has first been created 
from the coordinates of the edges and kinks and then refined to better match the actual atomic 
structure on the edges of the STM image. 

Atomic models of the GQDs were created for the TB calculations. For the smallest hexag- 
onal GQDs, this procedure was straight-forward as the hexagons can be readily identified by 
counting the carbon rings along each edge of the GQD. An example of the smallest observed 
GQD with its corresponding atomic model is shown in Supplementary Fig. 4(a). 

For the larger and more complex GQDs, a rough estimate of the structure was first created 
by measuring the coordinates of all the kinks and corners of the QGD from the STM image. 
An atomic model based on the coordinates was then created and compared to the atomically 
resolved image. As there might be an uncertainty in the STM piezo calibrations or some 
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residual thermal drift, the model was then further refined by comparing the number of 
carbon rings in the image and model along the edges of the GQD. These models correspond 
to the actual GQD geometry with an uncertainty of 1 carbon atom row. The kinks and 
other features on the edge are nevertheless produced with atomic precision. An example of 
this is shown in Supplementary Fig. 4(b). 

For solving the Klein-Gordon equation (which does not include atomic details), the 
shape and size were taken directly from the STM images for the larger QDs. The anal- 
ysis of the GQD size and shape was done using the free SPM analysis software Gwyddion 
(http://gwyddion.net/). The kinks on the edges of the larger flakes were also taken into 
account. The symmetric smaller QDs were modeled as perfect hexagons with sharp corners, 
where the edges run along the outermost atoms of the atomic model [hexagon surrounding 
the atomic model in Supplementary Fig. 4(a)]. 



dI/dV h MAPS FROM FIG. 2B WITH THE BACKGROUND 
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FIG. 5. Left: dln//dlnVb spectra measured in the middle of the GQD shown in the inset (blue 
solid line) and on top of the Ir(lll) substrate (red dashed line). This is the same GQD as shown 
in Fig. 2 of the main text. Right: Experimental dl /dV^ maps from Fig. 2b of the main text at 
the biases indicated in the figure showing the background, i.e. the surrounding Ir(lll) substrate. 
Set-point current 1 nA. 
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dI/dV h MAPS MEASURED AT POSITIVE BIASES 

topograph -290 mV 200 mV 




FIG. 6. dl/dVh maps measured at positive bias showing no contrast over the graphene QD. Map 
corresponding to the energy of the lS-state (bias voltage -0.29 V) is shown for comparison. Set- 
point current 0.5 nA. Topography of the GQD is shown on the top-right panel (bias voltage 0.1 V 
and set-point current 0.5 nA. 

CALCULATION OF THE LDOS FROM THE KLEIN-GORDON EQUATION 

The use of the Klein-Gordon equation to describe the quantum confined states in graphene 
QDs can be motivated starting from the Dirac equation [4] 

-iv F Ha -Vip = Eip (1) 

where a is the Pauli matrix and tp = (^) is wavefunction for the two sublattices a and b. 
Writing this out in component form gives 

-iv F h(d x - id y )^ b = Eip a 
-iv F H(d x + id y )ip a = Eip b 
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Separating the components gives the Klein-Gordon equation for both sublattices 



'a 



(3) 



v 2 h 2 \7 2 *jj b = E 2 



'A 



For an infinite graphene flake (i.e. ignoring the effect of boundary conditions), all the 
solutions of the Dirac equation are also solutions of the Klein-Gordon equation. 

The Klein-Gordon equation was solved numerically by the finite element method (FEM) 
using the commercial software Comsol Multiphysics (v. 3.5). Eigenstates of the Klein- 
Gordon equation were calculated up to 2 eV (using v F = 6.2 x 10 5 m/s). LDOS maps for a 
given energy E were created by weighing each eigenstate ipi by the energy difference of the 
eigenvalue Ei and E using a Gaussian distribution 



We checked the effect of the width a on the calculated LDOS maps: 60 meV was found to give 
the best match with the experimental results and was thus used for all of the calculations. 

TIGHT-BINDING CALCULATIONS 

We used a single-electron tight-binding (TB) model, suitable for describing the 7r-electrons 
of graphene. Hoppings up to third-nearest neighbors were included, with the original pa- 
rameters being t\ = —2.7 eV, t 2 = —0.20 eV and t 3 = —0.18 eV for first, second and 
third-nearest neighbors, respectively [9]. In order to match the experimental bias voltages, 
the parameters were scaled to t\ = —2.26 eV, t 2 = —0.168 eV and t 3 = —0.151 eV, keeping 
their ratios constant. In 2D graphene, these parameters correspond to a Fermi velocity of 
about 6.3 x 10 5 m/s. The LDOS was computed with a 60 meV Lorentzian broadening of 
the spectrum. In the figures, the LDOS from the TB model has been averaged over the sites 
in each hexagon of the graphene lattice. 

VISUAL COMPARISON OF MEASURED AND CALCULATED MAPS 

To find the correspondence between the measured and calculated states the LDOS maps 
were visually compared. This was done by calculating several LDOS maps in relatively 
small energy steps (6 meV) and comparing these to each of the measured dl/dV^ maps. 
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FIG. 7. LDOS maps calculated based on the KG equation. Map at the energy of 0.341 eV (green 
underline) is taken as the best match to the measured dl/dV^ map. The red markers denote the 
limits of error for the visual match. 

Since in most cases the calculated states change fairly little between such small energy 
steps, a range of possible matches was picked for each measured dl /dVh map. The energy 
in the middle of the range was used as the best match. The width of the range of the 
corresponding experimental and calculated LDOS maps was typically 50 meV. An example 
of this procedure is shown in Supplementary Fig. 7, where the best match is indicated by 
the green bar. 

Subsequently, the Fermi velocity has been scaled such that the slope of the plot of the 
theoretical energies vs. the experimental bias voltage is equal to unity. This gives the best- 
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fit Fermi velocity. In the case of the TB calculations (which has three hopping parameters), 
we scale all the hopping parameters such that their ratios remain constant. The intercept 
of the plot of the theoretical energies vs. the experimental bias voltages gives the doping of 
our graphene in the bulk limit. It is not a freely adjustable parameter but rather a result of 
the comparison between theory and experiment. 
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